This document presents the subset of the figures used for paper about monitoring.
library(dplyr)
library(tidyr)
library(ggplot2)
library(ggmap)
library(scatterpie)
library(rgdal)
library(multcompView)
library(car)
library(ggpmisc)
library(chisq.posthoc.test)
library(vcd)
Sampling data correspond to the data collected with kobo and previously cleaned with the script 1_preprocesamiento_datos_kobo.Rmd.
# load data
muestreo_tidy<-read.delim("../data/kobo/muestreo_dic2020_tidy.txt", header = TRUE)
parcelas_tidy<-read.delim("../data/kobo/parcelas_dic2020_tidy.txt", header = TRUE)
# pivot long parcelas data to have health data as a single variable
parcelas_long<-pivot_longer(parcelas_tidy,
cols = healthy:worm,
names_to = "tree_health_simplified",
values_to = "n_trees")
Data analyzed here correspond only to the trees that were approved during the validation by manually reviewing the photographs in kobotoolbox. Total of 1778 trees sampled, 1765 were approved in the validation.
muestreo_tidy<- filter(muestreo_tidy, X_validation_status=="validation_status_approved")
# write.csv(muestreo_tidy, file="../../../../Desktop/muestreo_tidy.csv")
Color palettes:
# Make a nice color pallete and legend order for all plots
my_cols=c("darkgreen",
"darkred",
"orangered1",
"cadetblue",
"tan",
"beige",
# "burlywood4",
"coral",
"aquamarine3",
"gray70",
"black")
desired_order=c("healthy",
"ozone",
"ozone_and_other",
"others_combined",
"drougth",
"fungi",
# "insect",
"worm",
"acid_rain",
"other",
"dead")
desired_names=c("healthy",
"ozone",
"ozone and other",
"others combined",
"drougth",
"fungi",
# "insect",
"worm",
"acid rain",
"other",
"dead")
# For ozone damage percentage
my_cols2<-c("darkgreen", "gold2", "chocolate1", "orangered", "red4", "darkorchid4")
desired_order_percentage<-c("0%","less than 10%", "10 to 40%", "40 to 50%", "50 to 70%", "more than 70%")
Multiplot fun:
# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols: Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
library(grid)
# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)
numPlots = length(plots)
# If layout is NULL, then use 'cols' to determine layout
if (is.null(layout)) {
# Make the panel
# ncol: Number of columns of plots
# nrow: Number of rows needed, calculated from # of cols
layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
ncol = cols, nrow = ceiling(numPlots/cols))
}
if (numPlots==1) {
print(plots[[1]])
} else {
# Set up the page
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
# Make each plot, in the correct location
for (i in 1:numPlots) {
# Get the i,j matrix positions of the regions that contain this subplot
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col))
}
}
}
Configure google api for maps:
# code adapted from https://rgraphgallery.blogspot.com/2013/04/rg-plot-pie-over-g0ogle-map.html
## configure google api
# You first need to register your api key in https://cloud.google.com/maps-platform/#get-started and follow instructions. The geocoding API is a free service, but you nevertheless need to associate a credit card with the account. Please note that the Google Maps API is not a free service. There is a free allowance of 40,000 calls to the geocoding API per month, and beyond that calls are $0.005 each.
# after you obtain your api, save it in /scripts/api_key.api (not shown in this repo por obvious reasons).
# if you get the following error when running get_map():
#"Error in aperm.default(map, c(2, 1, 3)) :
# invalid first argument, must be an array "
# check this troubleshooting: https://rgraphgallery.blogspot.com/2013/04/rg-plot-pie-over-google-map.html
## load and register api
api <- readLines("api_key.api")
register_google(key = api)
Map and monitoring figures presented in the paper:
# get cdmx shape
CDMX<-readOGR(dsn="../data/spatial", layer="CDMX")
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/veronicareyesgalindo/Documents/GitHub/monitoreo-oyameles/data/spatial", layer: "CDMX"
## with 1 features
## It has 8 fields
CDMX<-fortify(CDMX)
# get PNDL shape
PNDL<-readOGR(dsn="../data/spatial", layer="Desierto_Leones_Geo_ITRF08")
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum International_Terrestrial_Reference_Frame_2008 in
## Proj4 definition: +proj=longlat +ellps=GRS80 +no_defs
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/veronicareyesgalindo/Documents/GitHub/monitoreo-oyameles/data/spatial", layer: "Desierto_Leones_Geo_ITRF08"
## with 1 features
## It has 14 fields
PNDL<-fortify(PNDL)
# get background map
sat_map = get_map(location = c(lon = -99.133549, lat = 19.3), zoom = 10, maptype = 'terrain-background', source = "google")
## plot
p_a<-ggmap(sat_map) +
geom_polygon(data = CDMX,
aes(x = long, y = lat, group = group),
color="black", fill=NA, size=1.5) +
geom_polygon(data = PNDL,
aes(x = long, y = lat, group = group),
color="red", fill=NA, size=1.5) +
geom_point(aes(x=-98.95, y=19.6),
shape=0, stroke=2, size=5, color="black") +
geom_point(aes(x=-98.95, y=19.55),
shape=0, stroke=2, size=5, color="red") +
geom_text(aes(label="CDMX", x=-98.87, y=19.6),
color="Black", fontface="bold", size=5) +
geom_text(aes(label="PNDL", x=-98.87, y=19.55),
color="Black", fontface="bold", size=5) +
theme(text = element_text(size = 20))+
ggtitle("a)")
# get background map
sat_map = get_map(location = c(lon = -99.30, lat = 19.31), zoom = 13, maptype = 'satellite', source = "google")
## add towns names
towns<-data.frame(nombre=c("San Bartolo Ameyalco",
"Santa Rosa Xochiac",
"San Mateo Tlaltenango"),
long=c(-99.270, -99.29, -99.276),
lat=c(19.333, 19.325, 19.346))
## plot
p_b<-ggmap(sat_map) +
geom_polygon(data = PNDL,
aes(x = long, y = lat, group = group),
color="red", fill=NA, size=1.5) +
geom_point(data=towns, aes(x=long, y=lat), colour="red", size=1.5) +
geom_text(data=towns, aes(label=nombre, x=long, y=lat),
color="white", fontface="bold",
size=5, nudge_y=0.003) +
# add Cruz de Coloxtitla (CX), and Convento (Cn) landmarks
geom_text(aes(label="X", x=-99.3014, y=19.286068),
color="white", fontface="bold", size=4) +
geom_text(aes(label="C", x=-99.31, y=19.3133),
color="white", fontface="bold", size=4) +
theme(text = element_text(size = 20))+
ggtitle("b)")
## plot map
# get map
sat_map = get_map(location = c(lon = -99.3060, lat = 19.2909), zoom = 14, maptype = 'satellite', source = "google")
# plot sampled plots
p_c <- ggmap(sat_map)
p_c <- p_c + geom_point(data=parcelas_tidy,
aes(x=X_coordinates_longitude,
y=X_coordinates_latitude),
color="red") +
geom_text(data=parcelas_tidy,
aes(x=X_coordinates_longitude,
y=X_coordinates_latitude,
label=plot),
color="white",
check_overlap = TRUE,
hjust = 0, vjust=1, nudge_x = 0.0005,
size= 5) +
theme(text = element_text(size = 20))+
ggtitle("c)")
p_c